function [etrh, detrb, etri, etr2i] = hydro(X_GRID_SIZE, L, Z_GRID_SIZE, TIME_FINAL, TIME_STEPS, amp, fm, evalpoints, metricorder, ic, temp)

  TEMP_OFFSET = 1.0;
  TEMP_OFFSET = 150;
  TEMP_OFFSET = temp;

  global timer1 timer2 timer3 timer4 timer5;
  timer1 = 0;
  timer2 = 0;
  timer3 = 0;
  timer4 = 0;
  timer5 = 0;
  % package parameters
  params = h2params(TIME_FINAL, TIME_STEPS, X_GRID_SIZE, L, Z_GRID_SIZE, TEMP_OFFSET);

  params.timeOrder = 4;


  params.metricOn = false;
  params.metricMark = false;
  params.metricEvaluationPoints = evalpoints;
  params.metricLinearFnInterpolation = true;

  params.plotHydro = true;
  params.plotMetric = false;

  %if strcmp(ic, 'random')
  if findstr(ic, 'random') 
    params.hydroRandGradOn   = true;
  end
  params.hydroRandGradFreq = fm;

  if findstr(ic, 'wave')
    params.hydroPureWave     = true;
  end
  params.hydroPureWaveFreq = fm;
  params.hydroPureWaveAmp  = amp;

  params.hydroPeak         = false;
  params.hydroPeakAmp      = 0.1;
  params.hydroPeakRelWidth = 0.05;


  params.hydroVortices = false;
  params.hydroVorticesParams = [0.40 0.5 0.1 0.00005; 0.60 0.5 0.1 0.00005]

  params.hydroOrder = 2;

  params.hydroSaveOutput   = false;

  % summary
  disp(sprintf('Grid Size   : (%d, %d, %d)', params.Nx, params.Ny, params.Nz));
  disp(sprintf('Grid Length : %f', params.Lx));
  disp(sprintf('Cell Length : %f', params.Dx));
  disp(sprintf('Target Time : %f', params.Tfinal));
  disp(sprintf('Time steps  : %d',  params.Tsteps));
  disp(sprintf('Hydro Size  : %6.2fMB', 3*params.Nx.*params.Ny.*5.*TIME_STEPS.*64/8/1024/1024));
  disp(sprintf('Metric Size : %6.2fMB', length(params.metricEvaluationPoints).*params.Nx.*params.Ny.*params.Nz.*4.*4.*3.*64./8./1024./1024.*(params.timeOrder+1)));

  % prepare figures for plotting
  close all;

  % load color maps
  load('cm', 'mycmap');
  params.colormapOffset = 512;

  if (params.plotHydro || (params.plotMetric && params.metricOn))
    fig  = figure; 
    set(fig, 'Colormap', mycmap);
    set(fig, 'WindowStyle', 'docked');
  end
  refresh();
  pause(0.001);


  % initialize derivative kernels
  h2dinit(params);

  [tOut, uOut, udtOut, udtdtOut] = h2hydro(params);

  disp(sprintf('%.2f - derivatives', timer1));
  disp(sprintf('%.2f - m', timer2));
  disp(sprintf('%.2f - r', timer3));
  disp(sprintf('%.2f - inversion', timer4));
  disp(sprintf('%.2f - packaging', timer5));

  if (params.hydroSaveOutput)
    save(sprintf('%s.mat', datestr(clock, 'yy-mm-dd HHMM')), 'tOut', 'uOut', 'udtOut', 'udtdtOut', 'params', '-v7.3');
  end

  grerr = 0;
  if (params.metricOn == false)
    return
  end

  etrh   = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  detrb  = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  etri   = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  etr2i  = zeros(length(metricorder), length(evalpoints), params.Nx, params.Ny);
  

  for i=1:length(metricorder)
    params.metricOrder = metricorder(i);
    [g gt gtt] = h2metric(params, tOut, uOut, udtOut, udtdtOut);
    tic
    [etrh(i,:,:,:), detrb(i,:,:,:), etri(i,:,:,:), etr2i(i,:,:,:)] = h2gr(params, g, gt, gtt);
    disp(sprintf('h2gr in %.2f sec', toc))
  end

end


















